#include "triangle.hpp"
namespace chen {

bool Triangle::interact (const Ray & _ray, SurfInteract * si) const {
    Point3 A(a-_ray.o), B(b-_ray.o), C(c-_ray.o);
    Vector3  r_dir = _ray.d;
    int dim0 = 0;
    if(abs(r_dir[dim0]) < abs(r_dir[1])) dim0 = 1;
    if(abs(r_dir[dim0]) < abs(r_dir[2])) dim0 = 2; 
    if(dim0 == 0) {
        r_dir.permute(1, 2, 0);
        A.permute(1, 2, 0);
        B.permute(1, 2, 0);
        C.permute(1, 2, 0);
    } else if(dim0 == 1) {
        r_dir.permute(2, 0, 1);
        A.permute(2, 0, 1);
        B.permute(2, 0, 1);
        C.permute(2, 0, 1);
    }
    Float dxf = r_dir[0]/r_dir[2], dyf = r_dir[1]/r_dir[2];
    A.x -= dxf*A.z, A.y -= dyf*A.z;
    B.x -= dxf*B.z, B.y -= dyf*B.z;
    C.x -= dxf*C.z, C.y -= dyf*C.z;
    Float e0 = B.y*A.x-B.x*A.y;
    Float e1 = C.y*B.x-C.x*B.y;
    Float e2 = A.y*C.x-A.x*C.y;
#ifndef USE_DOUBLE_FLOAT
    if(e0 == 0 || e1 == 0 || e2 == 0) {
        double ed0 = (double)B.y*A.x-(double)B.x*A.y;
        double ed1 = (double)C.y*B.x-(double)C.x*B.y;
        double ed2 = (double)A.y*C.x-(double)A.x*C.y;
        if((ed0 < 0 || ed1 < 0 || ed2 < 0) && (ed0 > 0 || ed1 > 0 || ed2 > 0))
            return false;
        if(ed0 + ed1 + ed2 == 0) return false; // edge-on hit
    }
#endif
    if((e0 < 0 || e1 < 0 || e2 < 0) && (e0 > 0 || e1 > 0 || e2 > 0))
        return false;
    if(e0 + e1 + e2 == 0) return false; // edge-on hit
    Float dinv = 1/(e0+e1+e2);
    Float b0 = e0*dinv, b1 = e1*dinv, b2 = e2*dinv;
    Float tn = (C.z*b0+A.z*b1+B.z*b2)/r_dir[2];
    if(tn > _ray.t || tn < 0) return false;
    _ray.t = abs(tn);
    si->p = c*b0+a*b1+b*b2;
    Float xabss = abs(c.x*b0)+abs(a.x*b1)+abs(b.x*b2);
    Float yabss = abs(c.y*b0)+abs(a.y*b1)+abs(b.y*b2);
    Float zabss = abs(c.z*b0)+abs(a.z*b1)+abs(b.z*b2);
    si->p_err = Vector3(xabss, yabss, zabss)*gamma(7);
    si->ng = Normal(cross(b-a, c-a).normalized());
    si->ngy = (Normal)((b-a).normalized());
    si->flipped = dot((Vector3)si->ng, _ray.d) > 0;
    si->wh = _ray.d.normalized();
    si->hit = (const Object*)this;

    // Calculate uv coordinates and relative differentials
    Float numu0 = (b.y-a.y)*(c.x-a.x)-(b.x-a.x)*(c.y-a.y);
    Float numu1 = (b.z-a.z)*(c.x-a.x)-(b.x-a.x)*(c.z-a.z);
    Float numu2 = (b.z-a.z)*(c.y-a.y)-(b.y-a.y)*(c.z-a.z);
    Float u, v;
    if(abs(numu0) > abs(numu1) && abs(numu0) > abs(numu2)) {
        u = (si->p.y*(c.x-a.x)-(c.y-a.y)*si->p.x)/numu0;
        v = ((b.y-a.y)*si->p.x-(b.x-a.x)*si->p.y)/numu0;
    } else if(abs(numu1) > abs(numu0) && abs(numu1) > abs(numu2)) {
        u = (si->p.z*(c.x-a.x)-(c.z-a.z)*si->p.x)/numu1;
        v = ((b.z-a.z)*si->p.x-(b.x-a.x)*si->p.z)/numu1;
    } else {
        u = (si->p.z*(c.y-a.y)-(c.z-a.z)*si->p.y)/numu2;
        v = ((b.z-a.z)*si->p.y-(b.y-a.y)*si->p.z)/numu2;
    }
    Vector3 dpdu = b-a, dpdv = c-a;
    si->dpdu = dpdu;
    si->dpdv = dpdv;
    si->uv = {u, v};
    si->dndu = si->dndv = {0, 0, 0};
    return true; 
}
bool Triangle::mask(const Ray & _ray) const {
    Point3 A(a-_ray.o), B(b-_ray.o), C(c-_ray.o);
    Vector3  r_dir = _ray.d;
    int dim0 = 0;
    if(abs(r_dir[dim0]) < abs(r_dir[1])) dim0 = 1;
    if(abs(r_dir[dim0]) < abs(r_dir[2])) dim0 = 2; 
    if(dim0 == 0) {
        r_dir.permute(1, 2, 0);
        A.permute(1, 2, 0);
        B.permute(1, 2, 0);
        C.permute(1, 2, 0);
    } else if(dim0 == 1) {
        r_dir.permute(2, 0, 1);
        A.permute(2, 0, 1);
        B.permute(2, 0, 1);
        C.permute(2, 0, 1);
    }
    Float dxf = r_dir[0]/r_dir[2], dyf = r_dir[1]/r_dir[2];
    A.x -= dxf*A.z, A.y -= dyf*A.z;
    B.x -= dxf*B.z, B.y -= dyf*B.z;
    C.x -= dxf*C.z, C.y -= dyf*C.z;
    Float e0 = B.y*A.x-B.x*A.y;
    Float e1 = C.y*B.x-C.x*B.y;
    Float e2 = A.y*C.x-A.x*C.y;
#ifndef USE_DOUBLE_FLOAT
    if(e0 == 0 || e1 == 0 || e2 == 0) {
        double ed0 = (double)B.y*A.x-(double)B.x*A.y;
        double ed1 = (double)C.y*B.x-(double)C.x*B.y;
        double ed2 = (double)A.y*C.x-(double)A.x*C.y;
        if((ed0 < 0 || ed1 < 0 || ed2 < 0) && (ed0 > 0 || ed1 > 0 || ed2 > 0))
            return false;
        if(ed0 + ed1 + ed2 == 0) return false; // edge-on hit
    }
#endif
    if((e0 < 0 || e1 < 0 || e2 < 0) && (e0 > 0 || e1 > 0 || e2 > 0))
        return false;
    if(e0 + e1 + e2 == 0) return false; // edge-on hit
    Float dinv = 1/(e0+e1+e2);
    Float b0 = e0*dinv, b1 = e1*dinv, b2 = e2*dinv;
    Float tn = (C.z*b0+A.z*b1+B.z*b2)/r_dir[2];
    if(tn > _ray.t || tn < 0) return false;
    _ray.t = abs(tn);
    return true; 
}
void Triangle::uniformSample (const Point2 & _u, Interact * in, Point2 * uv, Float * pdf) const {
    Point2 u;
    if(_u.x +_u.y < 1) u = _u;
    else u = Point2(1-_u.x, 1-_u.y);
    in->p = a*u.x+(Vector3)(b*u.y)+(Vector3)(c*(1-u.x-u.y));
    in->ng = Normal(cross(b-a, c-a).normalized());
    in->p_err = (abs((Vector3)a)+abs((Vector3)b)+abs((Vector3)c))*gamma(5);
    in->flipped = false;
    in->hit = this;
    *in = this->w2o->reversed(*in);
    if(uv != nullptr) *uv = u;
    *pdf = (Float)1.0/area();
}
void Triangle::uniformSampleFrom (const Point2 & _u, const Interact & from, Interact * in, Point2* uv, Float * pdf) const {
    Point3 pfrom = (*(this->w2o))(from.p);
    Point2 u;
    if(_u.x +_u.y < 1) u = _u;
    else u = Point2(1-_u.x, 1-_u.y);
    in->p = a*u.x+(Vector3)(b*u.y)+(Vector3)(c*(1-u.x-u.y));
    in->ng = Normal(cross(b-a, c-a).normalized());
    in->p_err = (abs((Vector3)a)+abs((Vector3)b)+abs((Vector3)c))*gamma(5);
    in->flipped = dot((Vector3)in->ng, pfrom-in->p) < 0;
    in->hit = this;
    *in = this->w2o->reversed(*in);
    if(uv != nullptr) *uv = u;
    *pdf = 1/area();
}
}